Metropolis Hastings with EM MAP Start

library("spectralmc")

# get logprior,
# syntax: lwidth = prior.Gamma(shape, rate) means that lwidth has gamma prior with given shape and rate
logprior <- LogPrior(
  amp = prior.Positive(),
  lwidth = prior.Gamma(1.01, 0.01), #shape rate
  gwidth = prior.Gamma(1.01, 0.01), #shape rate
  pos = prior.Unif(0, 1000), # min, max
  a = prior.Norm(0, 0.00001), # mean, sd
  b = prior.Unif(-1000, 1000)
)

proposal <- Proposal(
  amp = prop.Norm(0.005), #sd
  lwidth = prop.Norm(0.008),
  gwidth = prop.Norm(0.008), # can also use prop.Unif(radius) for uniform proposal
  pos = prop.Norm(0.006),
  a = prop.Norm(0.00000005),
  b = prop.Norm(0.0002) #b = prop.Unif(0.02)
)

fit_kp <- 3 # number of peaks to fit
noise_sigma <- 0.001 # noise_sigma passed to likelihood
signal_model <- voigt.model # which signal to fit (could also be pseudovoigt or GoM)

# get synthetic data
data <- data.synthetic.stormyclouds(seed = 2021, kp = 3, noise = 0.001)
x <- data$x
y <- data$y

posterior <- Posterior(Loglikeli(signal_model, noise_sigma = noise_sigma), logprior)

start_params <- initial.fit.EM(x, y, fit_kp, seed=2021)

  |================================================================================| 100% elap:  0s eta:  0s

res <- chains.mhmc(x, y, start_params, posterior, proposal, iter = 20000, print_bar = FALSE)

  |================================================================================| 100% elap: 16s eta:  0s

Interactive plot to see fitting process (burn-in of chain): (Not very interesting if EM start is used.)

plt.fit.chain.interactive(x, y, res$samples, signal_model, step_width = 1000)

Double click on legend to zoom to trace.

plt.traceplot(res$samples, "pos")

plt.traceplot(res$samples, "amp")

plt.traceplot(res$samples, "lwidth")

plt.traceplot(res$samples, "gwidth")

plt.traceplot(res$samples, "a")

plt.traceplot(res$samples, "b")

plt.traceplot(res$samples, "b")

# can also do normal coda plots
library(coda)
samples_matrix <- utils.samples_to_matrix(res$samples) #can also filter variables with which=c("pos", "amp")
chain <- mcmc(samples_matrix)
summary(chain)

Iterations = 1:20000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 20000

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

               Mean        SD  Naive SE Time-series SE
pos[1]    2.117e+02 1.000e-02 7.074e-05      9.450e-04
pos[2]    5.365e+02 3.752e-02 2.653e-04      1.345e-02
pos[3]    8.666e+02 2.550e-02 1.803e-04      6.722e-03
gwidth[1] 2.780e+00 4.322e-02 3.056e-04      6.808e-03
gwidth[2] 7.517e+00 3.124e-01 2.209e-03      1.627e-01
gwidth[3] 3.975e+00 1.315e-01 9.297e-04      5.933e-02
lwidth[1] 3.740e-01 2.984e-02 2.110e-04      6.793e-03
lwidth[2] 5.862e+00 2.538e-01 1.795e-03      2.774e-01
lwidth[3] 3.445e+00 7.564e-02 5.349e-04      3.708e-02
amp[1]    3.850e+00 1.910e-02 1.351e-04      2.596e-03
amp[2]    4.630e+00 4.187e-02 2.960e-04      1.497e-02
amp[3]    2.363e+00 2.697e-02 1.907e-04      9.631e-03
a         1.911e-06 2.108e-07 1.491e-09      5.145e-08
b         1.938e+01 9.823e-05 6.946e-07      1.267e-05

2. Quantiles for each variable:

               2.5%       25%       50%       75%     97.5%
pos[1]    2.117e+02 2.117e+02 2.117e+02 2.117e+02 2.117e+02
pos[2]    5.364e+02 5.365e+02 5.365e+02 5.366e+02 5.366e+02
pos[3]    8.665e+02 8.665e+02 8.666e+02 8.666e+02 8.666e+02
gwidth[1] 2.748e+00 2.766e+00 2.775e+00 2.785e+00 2.818e+00
gwidth[2] 7.136e+00 7.251e+00 7.489e+00 7.660e+00 8.415e+00
gwidth[3] 3.838e+00 3.914e+00 3.954e+00 3.992e+00 4.466e+00
lwidth[1] 3.127e-01 3.631e-01 3.747e-01 3.907e-01 4.190e-01
lwidth[2] 5.468e+00 5.660e+00 5.850e+00 6.093e+00 6.262e+00
lwidth[3] 3.332e+00 3.380e+00 3.428e+00 3.496e+00 3.588e+00
amp[1]    3.829e+00 3.842e+00 3.848e+00 3.855e+00 3.874e+00
amp[2]    4.562e+00 4.602e+00 4.626e+00 4.655e+00 4.750e+00
amp[3]    2.334e+00 2.346e+00 2.359e+00 2.370e+00 2.454e+00
a         1.195e-06 1.827e-06 1.941e-06 2.035e-06 2.206e-06
b         1.938e+01 1.938e+01 1.938e+01 1.938e+01 1.938e+01

plot(chain)

autocorr.plot(chain)